Kinetic boundary conditions in the lattice Boltzmann method 

Santosh Ansumali 
ETH-Ziirich, Department of Materials, 
Institute of Polymers 
ETH-Zentrum, Sonneggstr. 3, 
ML J 19, CH-8092 Zurich, Switzerland 



Iliya V. Karlin 

ETH-Ziirich, Department of Materials, 
Institute of Polymers 
ETH-Zentrum, Sonneggstr. 3, 
ML J 19, CH-8092 Zurich, Switzerland 

Abstract 

Derivation of the lattice Boltzmann method from the continuous kinetic theory [X. He and L. 
S. Luo, Phys. Rev. E 55, R6333 (1997); X. Shan and X. He, Phys. Rev. Lett. 80, 65 (1998)] is 
extended in order to obtain boundary conditions for the method. For the model of a diffusively 
reflecting moving solid wall, the boundary condition for the discrete set of velocities is derived, and 
the error of the discretization is estimated. Numerical results are presented which demonstrate 
convergence to the hydrodynamic limit. In particular, the Knudsen layer in the Kramers' problem 
is reproduced correctly for small Knudsen numbers. 

PACS numbers: 05.70.Ln, 47.11. +j 
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I. INTRODUCTION 



In recent years, the lattice Boltzmann method (LBM) has emerged as an alternative 
tool for the computational fluid dynamics 0. Originally, the LBM was developed as a 
modification of the lattice gas model [Q] . Later derivations []3], |] revealed that the method is 
a special discretiszation of the continuous Boltzmann equation. The derivation of the LBM 
0] from the Boltzmann equation is essentially based on Grad's moment method || , together 
with the Gauss-Hermite quadrature in the velocity space. 

Another important issue was to retain positivity of discrete velocities populations in 
the bulk. Recently, a progress has been achieved in incorporating the if -theorem into the 
method || 0, and thus retaining positivity of the populations in the bulk. On the 
contrary, despite of several attempts || |10[ [TT], [12], [13], |L4|, [15], []1| a fully consistent theory of 
the boundary condition for the method is still lacking. It appears that the concerns about 
positivity of the population, and the connection with the continuous case, are somewhat 
ignored while introducing the boundary condition. The way the no-slip condition for the 
moving wall is incorporated in the method [[U], [TTJ, |12| is especially prone to danger of loss 
of positivity of the populations at the boundary. A clear understanding of the boundary 
condition becomes demanding for the case of moving boundary, complicated geometries, 
chemically reactive or porous walls. 

The theory of boundary conditions for the continuous Boltzmann equation is sufficiently 
well developed to incorporate the information about the structure and the chemical processes 
on the wall |T7[. The realization that the LBM is a special discretization of the Boltzmann 
equation allows to derive the boundary conditions for the LBM from continuous kinetic 
theory. In this work we demonstrate how this can be done in a systematic way. 

The outline of the paper is as follows: In section II we give a brief description of the LBM. 
In section [III] we briefly describe how boundary condition is formulated for the continuous 
kinetic theory. In section |V| we derive the boundary condition for the LBM and in section 
|V]we demonstrate some numerical simulation to validate the result. 
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II. OVERVIEW OF THE METHOD 



In the LBM setup, one considers populations of discrete velocities Cj, where i = 1, . . . , b, 
at discrete time t. It is convenient to introduce 6-dimensional population vectors /. 
In the isothermal case considered below, local hydrodynamic variables are given as, 

b 

p = J^/iM), 

i=l 

The basic equation to be solved is 

/i(r + a, t + 1) - fi(r, t) = -Pa[f(r, t)]Ai[f(r, t)}, (2) 

where (3 is a fixed parameter in the interval [0, 1] and is related to the viscosity. A scalar 
function of the population vector a is the nontrivial root of the nonlinear equation 

H(f) = H(f + aA[f}). (3) 

The function a ensures the discrete-time if -theorem. In the previous derivations || [| of the 
LBM from the Boltzmann equation, a quadratic form for the equilibrium distribution func- 
tion / eq , was obtained by evaluating the Taylor series expansion of the absolute Maxwellian 
equilibrium on the nodes of a properly selected quadrature. This was done to ensure that 
the Navier-Stokes equation is reproduced up to the order 0(M 2 ), where M is the Mach 
number. However, the disadvantage of expanding equilibrium distribution function is that 
the condition of monotonicity of the entropy production is not guaranteed. In order to avoid 
this problem, in the entropic formulation [7], ^J, the Boltzmann H function, rather than 
the equilibrium distribution, is evaluated at the nodes of the given quadrature, to get the 
discrete version of the if -function as 

tf = X>mg), (4) 

where wi denotes the weight associated with the corresponding quadrature node Cj. In the 
Appendix [A], the derivation of the if-function is presented. Afterwards, the collision term 
is constructed from the knowledge of the if -function (Eq.(f|)). The collision term A is 
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constructed in such a way that it satisfies a set of admissibility conditions needed to have a 
proper if -theorem and conservation laws (see Ref. || for details). 

The LBM model with the BGK collision form |4], [18|], can be considered as a limiting 
case of the entropic formulation. To obtain the lattice BGK equation, the function a in 
the Eq. (^) is set equals to 2, and for the collision term A BGK form is chosen. The 
equilibrium function used in the BGK form is obtained as the minimizer of the if -function 
(Eq.(||)) subjected to the hydrodynamic constrains (Eq. ([[])), evaluated up to the order M 2 
0. Derivation of the boundary conditions done in the subsequent section applies to both 
the forms of the LBM. 



III. BOUNDARY CONDITION FOR THE BOLTZMANN EQUATION 



Following Ref. [117), we briefly outline how boundary condition is formulated in the 
continuous kinetic theory. We shall restrict our discussion to the case where the mass flux 
through the wall is zero. For the present purpose, a wall dR is completely specified at any 
point (r G dR) by the knowledge of the inward unit normal n, the wall temperature T w 
and the wall velocity Uw- Hereafter, we shall denote the distribution function in a frame of 
reference moving with the wall velocity as g(£), with £ = c — U w . The distribution function 
reflected from the non-adsorbing wall can be written explicitly, if the scattering probability 
is known. In explicit form, 

|£ ■ n\g(£, t)= [ \£ • n\g(£, t)B (£' - £) (£ • n > 0), (5) 

j£'n<0 

where the non-negative function B (£ — > £) denotes the scattering probability from the 
direction £' to the direction £. If the wall is non-porous and non- adsorbing, the total 
probability for an impinging particle to be re-emitted is unity: 

/ B(£ , -*)de = l- (6) 

J£n>0 

Eq. (H) and Eq. @ ensure that the reflected distribution functions are positive and the 
normal flux through the wall is zero. A further restriction on the form of function B is 



dictated by the condition of detailed balance IT? 



[£' ■ n\g e ^', p w , 0, T W )B {£' - £) = |£ • n|^(-|, p w , 0, T W )B (-£ ^ -£') . (7) 
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A consequence of this property is that, if the impinging distributions are wall-Maxwellian, 
then the reflected distributions are also wall-Maxwellian. Thus, 



n|/ q (-£, Pw, 0, T w ) = / |£' • n\g e ^', p w , 0, T W )B (£' - £) dg. (8) 

J£'n<0 

This equation can also be understood as a weaker statement of the detailed balance condition 
TTfl . This form of the detailed balance is very attractive for our present purpose because of 
its integral nature, so that a discretization can be done in a natural way. 

In this paper, we only consider the diffusive boundary conditions because the steps as- 
sociated with the discretization are easier to appreciate due to the mathematical simplicity 
in this case. In this model of the wall it is assumed that the out-going stream has com- 
pletely lost its memory about the incoming stream. Thus, the scattering probability B is 
independent of the impinging directions, and is equals to 



B (f - t) = r " ^- Wr' 'n Z' \ j>/ = B (i) . (9) 







•n|£ ec l(-£,pw,0,T w ) 






n 


^q(^,pw,0,T w )^' 



Thus, the explicit expression for the reflected distribution function is 
Li ^ l£ ' n|o(£', t)d£ 

9& t) = 7 eq^ nT W^ eq (~£, Pw, 0, T w ), (£ • n > 0). (10) 

We need to transform this equation into the stationary co-ordinate system. As the 
equilibrium distribution depends only on the difference between the particle velocity and 
the local velocity, we have 

/M) = j ^ n< °' ' . T w .J eq (c,Pw^w,T w ), ((c-Uw) • n > 0). 

(11) 

In the next section, we will show how the discretization of the equation ([□]) can be 
performed. 

IV. DISCRETIZATION OF THE BOUNDARY CONDITION 

In the derivation of the lattice Boltzmann equation for the bulk various integrals, which 
are evaluated at the nodes of a Gauss-Hermite quadrature, || are of the form 

1=1 exp (-£ 2 )0(£)d£, (12) 

J£<LR D 
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where D is the spatial dimension. This form of the integral is well approximated by the 
Gauss-Hermite quadrature. However, the situation is different on the boundary because 
integrals appearing in Eq. (|11|) are over half-space. The choice of the quadrature in the 
bulk was based on the properties of integrals in the R D . If we would evaluate the integrals 
in Eq. flTT| ) using a quadrature defined in the half-space, this may introduce an undesirable 
mismatch of the nodes of the quadrature used on the boundary and that in the bulk. Thus, 
we here apply the quadrature used in the bulk even for the boundary nodes. Next, we shall 
estimate the extra error introduced by this procedure in comparison to the discretization 
error present in the bulk. 

The discrete distribution function used in the LBM is the projection of the continuous 
distribution function in a finite dimensional orthonormal Hermite basis 0]. The equilibrium 
also need to be projected in this basis to have correct conservation laws. This solution has a 
major drawback that the positivity of the distribution function is lost in the truncation. This 
problem is circumvented, if we evaluate the Boltzmann if -function, rather than its minimizer 
under the constrains of conservation of the hydrodynamic variables, for the discrete case. 
Indeed, the local equilibrium can also be written as 

f ec l(c, Pw , U w , T w ) = exp (a + (3 ■ c + 7c 2 ) , (13) 

where a, (3 and 7 are the Lagrange multipliers needed for the minimization of the Boltzmann 
if -function under the constrains of conservation of the hydrodynamic variables. These La- 
grange multipliers are calculated from the requirement that the moments of the equilibrium 
distribution / ec l are known hydrodynamic quantities. Now, once we have evaluated the 
projection of the Boltzmann if-function on a finite dimensional Hermite basis, we calculate 
the equilibrium from the knowledge of the discrete if-function. It turns out that the equi- 
librium corresponding to the discrete if-function also has the same functional form as Eq. 
(|l3l). Only difference is that the Lagrange multipliers has to be calculated from the discrete 
conservation laws. One example of explicit form of such equilibrium distribution function is 
given in Ref. 0. 

First projecting the distribution functions in the Hermite basis and then evaluating the 
integrals appearing in Eq. ([TT|) by quadrature, we have 

f(ci,t) = — ^— ~ ; -/ eq (ci,£/w,Pw), ((c i -U w )-n>0), (14) 

££.n<ol(£i ' n )\f q ( c v U w,Pw) 
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where 

fM=wip Atilp) (15) 

denotes the rescaled distribution function evaluated at nodes of the quadrature. This 
rescaled distribution function is the distribution function used in the LBM |3], M. The 
discrete equilibrium distribution function / ec l is the projection of the equilibrium distri- 
bution on the finite dimensional Hermite basis and calculated by the procedure discussed 
above. Before estimating error associated with this formula, a few remark about the preced- 
ing equation is in order. First, in the isothermal case the wall temperature is a redundant 
quantity and is dropped from the argument of equilibrium distribution. To get a boundary 
condition for the lattice BGK equation, the true discrete equilibrium appearing in the Eq. 
([U]) can be replaced by the equilibrium used in the BGK model fll8| . This substitution 



is justified because up to order 0(M 2 ) the true equilibrium can be replaced by the BGK 
equilibrium. However, positivity of the reflected distributions may be lost in this truncation 
in the same way it happens in the bulk for lattice BGK model. A similar expression for the 
boundary conditions was earlier postulated by Gatingnol in the context of discrete velocity 
models of the kinetic theory PT . 



In order to estimate the extra error introduced on the boundary in comparison to the 
bulk, we write the ratio of two integrals appearing in the eq.flTTD as 

^■ n< ol^-n|/(^t)^ 

4n<0^'- n l/ eq ( C ''Pw,^W,T w X 

f,, n<0 \Z'-n\f^(c>,t)dc> (16) 
4„<o I*' • n l / eq ( c '> Pw, U w , T w )dg ' 

where f ne< i = f — f e< ^. As discussed above, the evaluation of the integral appearing in the 
denominator is straightforward in the sense that the order of accuracy of this evaluation is 
same as that of moments evaluation in the bulk. In order to evaluate the integral appearing in 
the numerator, we perform Hermite expansion of the non-equilibrium part of the distribution 
function around the zero velocity equilibrium. The result is 



N W 

f^(c } t) = f^(c, p w , 0, T w ) (c) • (17) 

i=o l ' 
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The first two expansion coefficients of the non-equilibrium part of the population are (a^ = 
0,af } = 0). In case of the isothermal hydrodynamics, only non-zero Hermite coefficient 

(2) 

needs to be kept is a\j . This is a symmetric tensor and is independent of the particle 
velocity c. After using the symmetry of the second order Hermite-polynomials, 

2 4 n<0 ^-n|/ e( l(c',pw,t/w,r w )rfc' • 1 j 

This expression can be evaluated using the Gauss-Hermite quadrature. The result is 

Et>.n<o\£i-*\f eq (c',PW,U w y 



This expression gives an estimate of the order of the accuracy of the eq. (|I4|) . In evaluation 
of the moments ( up to the second order moment) of the distribution function, no extra error 
is introduced as compared to the bulk. This happened because first odd order Hermite coef- 
ficient appearing in the expansion is zero. Due to the expansion around global equilibrium, 
used in the derivation, the boundary condition is valid only up to the order 0(M 2 ). 

Now, for purely diffusive scattering, we have a closed form expression for the reflected 
populations with the same order of accuracy as the bulk node. However, we have said 
nothing about the grazing directions. Unlike continuous kinetic theory, here we need to 
specify the conditions in the grazing directions. Only information we have about the grazing 
populations is their positivity. A simple way to fixed the grazing population is to let them 
evolve according to the lattice-Boltzmann equation like nodes in the bulk region. This 
condition is implemented in the simulations presented in the next section. 



V. NUMERICAL TESTS 

The boundary condition derived in the previous section (Eq. fll4|)), retains one important 
feature of original Boltzmann equation, the Knudsen number dependent slip at the wall. To 
show this, we have performed a numerical simulation of the Kramers' problem [|17[]. This 
is one of the few problem where solution of the continuous Boltzmann equation is known 
analytically. This problem is a limiting case of the plane Couette flow, where one of the plate 
is moved to infinity, while keeping a fixed shear rate. We compare the analytical solution 
for the slip-velocity at the wall calculated for the linearized BGK collision model with the 
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numerical solution in the Fig. [l|. We have performed the numerical computation for the 
D2Q9 lattice with the entropic formulation of the LBM [7], |8|] with the expression of the H 
function given by Eq. (|^5|). The agreement between the two result for Knudsen number 
going to zero is very good. This is indeed an important result as it shows that with the 
proper implementation of the boundary condition, the solution of the LBM converge to the 
hydrodynamic limit (Knudsen number going to zero) in the same way as the Boltzmann 
equation. 

By simulating the Kramer's problem, we have shown that the present boundary condition 
can be used for stationary wall. To validate the boundary condition for the moving wall, 
we have performed simulation of the lid-driven cavity flow. The plot of stream-function is 
given in Fig.|2] for Reynolds number Re = 1000. The location of the primary and secondary 
vortex and the magnitude of the stream-functions agrees well with the previous simulations 

23. 



Once we have shown that the diffusive boundary condition used for the continuous Boltz- 
mann equation can be reformulated for the discrete case, the question arises that, can this 
procedure be applied for a more sophisticated scattering kernels used in the continuous ki- 
netic theory (see for examples Ref. |17|]). The answer is in affirmative for any condition 
written in the integral form, while in general it cannot be done for a point-wise condition 
like purely specular reflection. For example, a very general form of scattering probability, 
written in the integral form and can be easily modified for nonzero mass flux, is given in the 
Ref. 0] (see Eq.(6.26) of the Chapter 3). This form of the scattering probability can be 
discretized using the present method. 

It is instructive to compare the 'bounce-back' condition used in the literature with the 
present boundary condition. It can be seen easily that the present boundary condition re- 
duces to the bounce-back condition for the three velocity model used in the one-dimensional 
case. However, there is no correspondence between the two condition in the higher dimen- 
sions. 

The present boundary condition retains the positivity at the boundary nodes. This is a 
major advantage in comparison to other proposed boundary conditions for the purpose of 
the numerical stability. The Knudsen number dependent wall slip is a manifestation of the 
kinetic nature of the lattice Boltzmann equation. This nature of the scheme can be a burden 
if one is interested in solving the macroscopic creeping flow problems with very small grid 
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size. This will put some restriction on the simulation of creeping flow in very small grids 
(lattice Knudsen number = v/Lc s ). However, the restriction is not severe because of the fact 
that we still have the freedom to choose velocity very small to attain zero Reynolds number 
situation. In fact, the same condition is required for the validity of the LBM simulation of 
the hydrodynamics in the bulk. To conclude, we have proposed boundary conditions based 
on the kinetic theory considerations. A systematic way of dealing with the conditions at the 
boundary is developed for the lattice-Boltzmann method. The present work opens the way 
to the future development for the cases of reactive, porous or adsorbing walls. 
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APPENDIX A: DERIVATION OF THE if FUNCTION 



The Boltzmann if -function is 

H = J f log fd£ (Al) 

This expression written in terms of the logarithm of the distribution function (/i = log /) is 

H = J AtexpO) d$, (A2) 

We have chosen to work with the variable /i because the projection of it on to the Hermite 
basis preserves positivity of the distribution function /. The expansion of the function \i is 

Im = A<°>W<°> (c) + A^H® (c) + A%H% (c) (A3) 

This expression can also be written as, 

H = A<°>«<°> (c) + (c) + B§H% (c) - ^ (A4) 

The expansion coefficients A is calculated by the requirement that the moments of exp (//) 
are hydrodynamic variables. The expansion used here is a slightly different form of the 



Grad's moment expansion [Tj| and is known as the maximum entropy approximation [22 



23, 24]. Now, the Boltzmann if function is in a integral form suited for the evaluation in 
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the Gauss-Hermite quadrature (see Eq. (|12|)). This evaluation gives the discrete form of 
the H function as 

b 



(A5) 



i=l 



where, 



(A6) 



where To is the reference temperature. In 2-dimension, the nodes of the quadrature and the 
corresponding weights are 



and 



{0,0} 



7r(i-l) 
2 



2 



cv^(cos(« , sin 



if i = 



if i = 1,2,3,4 



if i = 5, 6, 7, 8, 



(A7) 



Wi 



± 

9 

1 
9 

J_ 

^ 36 



if i = 

if i= 1,2,3,4 

if i = 5,6,7,8. 



(AS 



Here the magnitude of the discrete velocity c is related to the reference temperature by 
the relation c = a/ (3/cbTq). With this, the entropy expression derived here coincide with 
the expression derived in Ref. B, by a different argument. 
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FIG. 1: Relative slip observed at the wall in the simulation of the Kramers' problem for shear 
rate a = 0.001, box length L = 32, Voo = a x L = 0.032. All the quantities are given in the 
dimensionless lattice unit. 
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FIG. 2: Stream-function for Re = 1000 in a simulation of lid driven cavity flow. Parameters used 
are: grid size 320 x 320, and lid velocity V = 0.075. All the quantities are given in the dimensionless 
lattice unit. 
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